home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2001 May / SGI IRIX Base Documentation 2001 May.iso / usr / share / catman / p_man / cat3 / complib / dpbsvx.z / dpbsvx
Encoding:
Text File  |  1998-10-30  |  11.4 KB  |  331 lines
  1.  
  2.  
  3.  
  4. DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DPBSVX - use the Cholesky factorization A = U**T*U or A = L*L**T to
  10.      compute the solution to a real system of linear equations  A * X = B,
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE DPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, EQUED,
  14.                         S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK,
  15.                         INFO )
  16.  
  17.          CHARACTER      EQUED, FACT, UPLO
  18.  
  19.          INTEGER        INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
  20.  
  21.          DOUBLE         PRECISION RCOND
  22.  
  23.          INTEGER        IWORK( * )
  24.  
  25.          DOUBLE         PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
  26.                         BERR( * ), FERR( * ), S( * ), WORK( * ), X( LDX, * )
  27.  
  28. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  29.      DPBSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
  30.      compute the solution to a real system of linear equations
  31.         A * X = B, where A is an N-by-N symmetric positive definite band
  32.      matrix and X and B are N-by-NRHS matrices.
  33.  
  34.      Error bounds on the solution and a condition estimate are also provided.
  35.  
  36.  
  37. DDDDEEEESSSSCCCCRRRRIIIIPPPPTTTTIIIIOOOONNNN
  38.      The following steps are performed:
  39.  
  40.      1. If FACT = 'E', real scaling factors are computed to equilibrate
  41.         the system:
  42.            diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
  43.         Whether or not the system will be equilibrated depends on the
  44.         scaling of the matrix A, but if equilibration is used, A is
  45.         overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
  46.  
  47.      2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
  48.         factor the matrix A (after equilibration if FACT = 'E') as
  49.            A = U**T * U,  if UPLO = 'U', or
  50.            A = L * L**T,  if UPLO = 'L',
  51.         where U is an upper triangular band matrix, and L is a lower
  52.         triangular band matrix.
  53.  
  54.      3. The factored form of A is used to estimate the condition number
  55.         of the matrix A.  If the reciprocal of the condition number is
  56.         less than machine precision, steps 4-6 are skipped.
  57.  
  58.      4. The system of equations is solved for X using the factored form
  59.         of A.
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))
  71.  
  72.  
  73.  
  74.      5. Iterative refinement is applied to improve the computed solution
  75.         matrix and calculate error bounds and backward error estimates
  76.         for it.
  77.  
  78.      6. If equilibration was used, the matrix X is premultiplied by
  79.         diag(S) so that it solves the original system before
  80.         equilibration.
  81.  
  82.  
  83. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  84.      FACT    (input) CHARACTER*1
  85.              Specifies whether or not the factored form of the matrix A is
  86.              supplied on entry, and if not, whether the matrix A should be
  87.              equilibrated before it is factored.  = 'F':  On entry, AFB
  88.              contains the factored form of A.  If EQUED = 'Y', the matrix A
  89.              has been equilibrated with scaling factors given by S.  AB and
  90.              AFB will not be modified.  = 'N':  The matrix A will be copied to
  91.              AFB and factored.
  92.              = 'E':  The matrix A will be equilibrated if necessary, then
  93.              copied to AFB and factored.
  94.  
  95.      UPLO    (input) CHARACTER*1
  96.              = 'U':  Upper triangle of A is stored;
  97.              = 'L':  Lower triangle of A is stored.
  98.  
  99.      N       (input) INTEGER
  100.              The number of linear equations, i.e., the order of the matrix A.
  101.              N >= 0.
  102.  
  103.      KD      (input) INTEGER
  104.              The number of superdiagonals of the matrix A if UPLO = 'U', or
  105.              the number of subdiagonals if UPLO = 'L'.  KD >= 0.
  106.  
  107.      NRHS    (input) INTEGER
  108.              The number of right-hand sides, i.e., the number of columns of
  109.              the matrices B and X.  NRHS >= 0.
  110.  
  111.      AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
  112.              On entry, the upper or lower triangle of the symmetric band
  113.              matrix A, stored in the first KD+1 rows of the array, except if
  114.              FACT = 'F' and EQUED = 'Y', then A must contain the equilibrated
  115.              matrix diag(S)*A*diag(S).  The j-th column of A is stored in the
  116.              j-th column of the array AB as follows:  if UPLO = 'U',
  117.              AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j; if UPLO = 'L',
  118.              AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).  See below for
  119.              further details.
  120.  
  121.              On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
  122.              diag(S)*A*diag(S).
  123.  
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))
  137.  
  138.  
  139.  
  140.      LDAB    (input) INTEGER
  141.              The leading dimension of the array A.  LDAB >= KD+1.
  142.  
  143.      AFB     (input or output) DOUBLE PRECISION array, dimension (LDAFB,N)
  144.              If FACT = 'F', then AFB is an input argument and on entry
  145.              contains the triangular factor U or L from the Cholesky
  146.              factorization A = U**T*U or A = L*L**T of the band matrix A, in
  147.              the same storage format as A (see AB).  If EQUED = 'Y', then AFB
  148.              is the factored form of the equilibrated matrix A.
  149.  
  150.              If FACT = 'N', then AFB is an output argument and on exit returns
  151.              the triangular factor U or L from the Cholesky factorization A =
  152.              U**T*U or A = L*L**T.
  153.  
  154.              If FACT = 'E', then AFB is an output argument and on exit returns
  155.              the triangular factor U or L from the Cholesky factorization A =
  156.              U**T*U or A = L*L**T of the equilibrated matrix A (see the
  157.              description of A for the form of the equilibrated matrix).
  158.  
  159.      LDAFB   (input) INTEGER
  160.              The leading dimension of the array AFB.  LDAFB >= KD+1.
  161.  
  162.      EQUED   (input or output) CHARACTER*1
  163.              Specifies the form of equilibration that was done.  = 'N':  No
  164.              equilibration (always true if FACT = 'N').
  165.              = 'Y':  Equilibration was done, i.e., A has been replaced by
  166.              diag(S) * A * diag(S).  EQUED is an input argument if FACT = 'F';
  167.              otherwise, it is an output argument.
  168.  
  169.      S       (input or output) DOUBLE PRECISION array, dimension (N)
  170.              The scale factors for A; not accessed if EQUED = 'N'.  S is an
  171.              input argument if FACT = 'F'; otherwise, S is an output argument.
  172.              If FACT = 'F' and EQUED = 'Y', each element of S must be
  173.              positive.
  174.  
  175.      B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
  176.              On entry, the N-by-NRHS right hand side matrix B.  On exit, if
  177.              EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwritten
  178.              by diag(S) * B.
  179.  
  180.      LDB     (input) INTEGER
  181.              The leading dimension of the array B.  LDB >= max(1,N).
  182.  
  183.      X       (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
  184.              If INFO = 0, the N-by-NRHS solution matrix X to the original
  185.              system of equations.  Note that if EQUED = 'Y', A and B are
  186.              modified on exit, and the solution to the equilibrated system is
  187.              inv(diag(S))*X.
  188.  
  189.      LDX     (input) INTEGER
  190.              The leading dimension of the array X.  LDX >= max(1,N).
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))
  203.  
  204.  
  205.  
  206.      RCOND   (output) DOUBLE PRECISION
  207.              The estimate of the reciprocal condition number of the matrix A
  208.              after equilibration (if done).  If RCOND is less than the machine
  209.              precision (in particular, if RCOND = 0), the matrix is singular
  210.              to working precision.  This condition is indicated by a return
  211.              code of INFO > 0, and the solution and error bounds are not
  212.              computed.
  213.  
  214.      FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  215.              The estimated forward error bound for each solution vector X(j)
  216.              (the j-th column of the solution matrix X).  If XTRUE is the true
  217.              solution corresponding to X(j), FERR(j) is an estimated upper
  218.              bound for the magnitude of the largest element in (X(j) - XTRUE)
  219.              divided by the magnitude of the largest element in X(j).  The
  220.              estimate is as reliable as the estimate for RCOND, and is almost
  221.              always a slight overestimate of the true error.
  222.  
  223.      BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  224.              The componentwise relative backward error of each solution vector
  225.              X(j) (i.e., the smallest relative change in any element of A or B
  226.              that makes X(j) an exact solution).
  227.  
  228.      WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
  229.  
  230.      IWORK   (workspace) INTEGER array, dimension (N)
  231.  
  232.      INFO    (output) INTEGER
  233.              = 0:  successful exit
  234.              < 0:  if INFO = -i, the i-th argument had an illegal value
  235.              > 0:  if INFO = i, and i is
  236.              <= N: the leading minor of order i of A is not positive definite,
  237.              so the factorization could not be completed, and the solution has
  238.              not been computed.  = N+1: RCOND is less than machine precision.
  239.              The factorization has been completed, but the matrix is singular
  240.              to working precision, and the solution and error bounds have not
  241.              been computed.
  242.  
  243. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  244.      The band storage scheme is illustrated by the following example, when N =
  245.      6, KD = 2, and UPLO = 'U':
  246.  
  247.      Two-dimensional storage of the symmetric matrix A:
  248.  
  249.         a11  a12  a13
  250.              a22  a23  a24
  251.                   a33  a34  a35
  252.                        a44  a45  a46
  253.                             a55  a56
  254.         (aij=conjg(aji))         a66
  255.  
  256.      Band storage of the upper triangle of A:
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPBBBBSSSSVVVVXXXX((((3333FFFF))))
  269.  
  270.  
  271.  
  272.          *    *   a13  a24  a35  a46
  273.          *   a12  a23  a34  a45  a56
  274.         a11  a22  a33  a44  a55  a66
  275.  
  276.      Similarly, if UPLO = 'L' the format of A is as follows:
  277.  
  278.         a11  a22  a33  a44  a55  a66
  279.         a21  a32  a43  a54  a65   *
  280.         a31  a42  a53  a64   *    *
  281.  
  282.      Array elements marked * are not used by the routine.
  283.  
  284.  
  285.  
  286.  
  287.  
  288.  
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.